MOTION COMPENSATION IN BIQMAGNETIC 



MEASUREMENT 

Technical Field ; 
5 [0001] The invention relates to the measurement of biomagnetic 
signals. The invention has application to the magnetic imaging of 
biological structures. Some embodiments of the invention are used to 
compensate for head motion in magnetoencephalography. 

10 Background 

[0002] Magnetoencephalography ("MEG") involves detecting 
magnetic fields produced within a subject's brain. Information about 
such biomagnetic fields is most useful when it can be associated with 
particular structures within the subject's brain to a spatial resolution of a 

1 5 few millimeters or better. 

[0003] A typical MEG system comprises a helmet which carries a 
large number of magnetic detectors. The subject's head is placed inside 
the helmet. Magnetic signals originating from within the subject's head 
20 are detected at each of the magnetic detectors over a data acquisition 
period. The data acquisition period may, for example, be a few minutes. 

[0004] Since biomagnetic signals are typically measured over 
significant time spans it is unreasonable to expect that a subject will be 
25 able to hold his or her head completely still throughout the measurement. 
Motions of the subject's head during the data acquisition period can 
interfere with the ability to associate particular magnetic signals with 
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specific structures within the subject's brain. It is not practical to 
completely immobilize the subject's head. 



[0005] One can localize the subject's head relative to the measured 
5 magnetic fields if the position of the subject's head relative to the 

magnetic detectors used to detect the magnetic fields is known. One way 
to localize a subject's head is to attach small coils at three or more known 
locations on the subject's head. The coils create fluctuating magnetic 
fields when altemating electrical currents are passed through them. 
10 Magnetic detectors are used to detect the coils' magnetic fields. The 
location of the coils, and thus the location of the subject's head, can be 
determined from the detected magnetic fields of the coils. 

[0006] There is a need for biomagnetic measurement systems which 
15 can compensate for motion of the part of a subject being studied. 

Summary of the Invention 

[0007] One aspect of this invention provides a method for magnetic 
imaging of an object. The method includes monitoring a magnetic field 

20 of sources in the object at a plurality of magnetic detectors to obtain a 
corresponding plurality of sensor outputs. A position of the object is 
monitored while monitoring the magnetic field of the sources. The 
magnetic field of the sources in the object is modeled as a gradient of a 
scalar potential. The scalar potential includes a sum of spherical 

25 harmonic functions each multiplied by a corresponding coefficient. The 
method includes compensating for changes in the position of the object 
by applying a transformation to the plurality of sensor outputs. The 
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transformation includes, at least in part, a spherical harmonic translation 
transformation. 



[0008] Further aspects of the invention and features of specific 
5 embodiments of the invention are described below. 

Brief Description of the Drawings 

[0009] In drawings which illustrate non-limiting embodiments of 
the invention, 

10 Figure 1 is a schematic diagram of a MEG system which is used to 

measure biomagnetic fields according to one embodiment of the 
invention; and. 

Figure 2 is a flow chart which illustrates a method for magnetic 
imaging of an object according to another embodiment of the invention. 

15 

Description 

[001 0] Throughout the following description, specific details are 
set forth in order to provide a more thorough understanding of the 
invention. However, the invention may be practiced without these 
20 particulars. In other instances, well known elements have not been 
shown or described in detail to avoid unnecessarily obscuring the 
invention. Accordingly, the specification and drawings are to be 
regarded in an illustrative, rather than a restrictive, sense. 

25 [001 1] Figure 1 shows a MEG system 20 according to one 

embodiment of the invention. System 20 includes a number of magnetic 
detectors 24 which are located to detect magnetic fields originating 
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within the head of a subject S. While Figure 1 represents a detector array 
by only a few detectors 24, a MEG system typically has a few hundred 
detectors 24. 

5 [00121 Detectors 24 may be SQUID detectors. Detectors 24 may 
comprise magnetic gradiometers and/or magnetometers. Signals from 
detectors 24 are appropriately conditioned by way of suitable signal 
conditioning circuitry 26. The processed signals are provided to a signal 
processing system 28. Signal processing system 28 generates output at an 
1 0 output device 30. The output may comprise, for example, a graphical 
display, a data file containing MEG data, or the like. 

[0013] MEG system 20 includes a head localization system 21 
capable of monitoring the position and orientation of the subject's head 

15 during the acquisition of MEG data. In the illustrated embodiment, head 
localization system 21 comprises a system as described in the co-pending 
and commonly owned patent application entitled METHOD AND 
APPARATUS FOR LOCALIZING BIOMAGNETIC SIGNALS being 
filed simultaneously herewith, which is hereby incorporated herein by 

20 reference. Any other suitable localization system, now known or 

developed in future which provides information about the position, or 
changes in position of a subject's head or other body part could be used 
as localization system 21. 

25 [0014] In the illustrated embodiment, system 21 includes coils 32, 
33 and 34. Coils 32, 33 and 34 are attached to the subject's head and are 
respectively driven with signals of frequencies //,/^ and by a controller 



36 which, in the illustrated embodiment includes oscillators 37, 38, and 
39. Driving signals for coils 32, 33 and 34 may be obtained in any 
suitable manner. 

5 [0015] System 21 monitors the locations of coils 32, 33 and 34. 
System 21 determines the position and orientation of the subject's head 
at various times based upon the locations of coils 32, 33 and 34. This 
may be done, for example, using methods described in the above-noted 
co-pending application. 

10 

[0016] During a data acquisition period, MEG system 20 monitors 
the outputs of magnetic detectors 24 and the position of the subject's 
head. This yields a sequence of samples of the outputs from each of 
magnetic detectors 24. System 20 generates a map of magnetic sources 
15 within the subject's head based upon the samples. System 20 

compensates for motion of the subject's head based upon head position 
information provided by localization system 21. 

[0017] Figure 2 is a flow chart illustrating a method 200 according 
20 to one embodiment of the invention. Method 200 may be used for 

magnetic imaging of an object having sources which produce magnetic 
fields, such as the head of subject S of Figure 1, In a preferred 
embodiment, method 200 is carried out by MEG system 20 of Figure 1. 
After method 200 begins at block 202, it obtains a plurality of outputs 
25 from detectors 24 at block 204. At block 206, method 200 obtains the 
object's position from localization system 21. At block 208, the 
magnetic field of the object is modelled as a gradient of a scalar potential 



10 
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comprising a sum of spherical harmonic functions each multiplied by a 
corresponding coefficient, as described further below. At block 210, 
method 200 compares the object's position obtained in block 206 with a 
reference position, and compensates for motion of the object by 
correcting the magnetic field modelled in block 208 to obtain a plurality 
of motion corrected sensor outputs. The compensation carried out in 
block 210 preferably comprises applying a transformation to the 
spherical harmonic functions, as described further below. The 
transformation includes a spherical harmonic translation transformation. 



[0018] If one chooses a volume which includes detectors 24 but 
excludes any magnetic sources, the magnetic field detected by detectors 
24 can be expressed as the gradient of a magnetic scalar potential Y. 
Within the volume, in which there are no magnetic sources, Y satisfies 
1 5 Laplace's equation: 

= 0 (1) 
Y can be expressed as a sum of spherical harmonic functions, for 
example, as follows: 

^ ^ /=l ;w=-//j=0 \^0y /=1 /«=-/«=0 ^ ^ 



20 where /, m and n are integer indices; a,^„ are coefficients; is a 
normalizing radius; and Y/" are given by: 

P^'"{cos0)cos{m(p) l>m>0 



(3) 



where is the Schmidt semi-normalized associated Legendre function 
given by: 



and PP is the standard associated Legendre function as defined, for 

example, in Abramowitz and Stegun, Handbook of Mathematical 
Functions, National Bureau of Standards, 1964. This formulation has the 
advantage that it does not require calculations involving complex number 
types and the normalization is the same for different terms. Spherical 
harmonic terms with n=0 are called intemal terms, while terms with n=\ 
are called extemal terms. Sources inside the subject generate intemal 
terms in the spherical harmonic expansion, while sources outside the 
volume containing the detectors generate extemal terms in the spherical 
harmonic expansion. 

[001 9] The spherical harmonic functions could be normalized in 
any suitable manner. The spherical harmonic functions may be 
expressed in other mathematically equivalent ways. For example, the 
spherical harmonic functions could be expressed in a format which uses 
complex variables as follows: 




(5) 
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[0020] Where there are Af magnetic detectors 24, the signals from 
any one of the magnetic detectors may be represented as where / is an 
integer index corresponding to the magnetic detector with 1 <i<.M. 

5 [0021] The method uses the values to determine values for a 
number N of the coefficients which yield a magnetic scalar potential 
which would produce the observed magnetic fields. A'^is larger than M 
and is chosen to be sufficiently large that the finite series: 

Yu^J^^) (6) 

10 is a good approximation to the true scalar potential of the magnetic field 
being modelled. In equation (6) spherical harmonic functions. 

A^is not so large that the problem of determining values for the 
coefficients becomes impractical. 

1 5 [0022] The N coefficients are related to the M sensor outputs B^ 
by an Mx TV matrix as follows: 

Bn^lL^m^^ « ^=^^ (7) 

The matrix L (the "forward solution") is determined by the parameters 
of magnetic detectors 24 (e.g. their gain, coil position and size, and 
20 gradiometer order). The elements of Z may be computed in advance 

based upon the known geometry and properties of magnetic detectors 24. 
A minimum norm algorithm may be used to identify a set of coefficients 
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such that equation (7) is satisfied and a normaHzation function is 
minimized. 



[0023] Various volume integrals of gradients of Y can be used as 
5 the normaHzation function for the scalar potential T. The standard 
magnetic energy is one such volume integral. The magnetic energy can 
be expressed as: 




(8) 



10 An alternate normalizing function is given by: 




(9) 



E2 is a better choice for use as the normalization function where first- 
order gradiometers are used to detect the magnetic field. If magnetic 
detectors 24 include a sufficient number of higher-order gradiometers, 

1 5 then total second or higher-order derivatives could be used as the 

normalizing functions. A linear combination of J?y, E2 and/or other total 
higher-order derivatives could also be used as the normalization 
function. In an example embodiment, magnetic detectors 24 comprise 
first order gradiometers and the normalization function E2 defined by 

20 equation (9) is used as the normalization function. 
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[0024] When ¥ is represented by a spherical harmonic expansion 
then Equation (8) may be expressed as a matrix equation involving the 

A^xA'^ energy matrix K as follows: 

^2 =Z^P^P.^. 



(10) 



10 



where p and q each represent a particular set of values {l,m,n} and the 
summation in Equation (10) is performed over those TV sets of values 
corresponding to the spherical harmonic functions to be included in 
the model. Where the volume of integration is a spherical shell having 
inner radius rj and outer radius z-^? the elements of K are given by two 
internal basis functions and two external basis functions. The two 
intemal basis functions are given by: 



^.,o.Mo = 4^/ + 1)(/ + 2) 



2/+3 



2/+3 



(11) 



The two external basis functions are given by: 



(12) 



1 5 where is a normahzing radius. The intemal functions («=0) and 

external functions («=1) are orthogonal and so A'„/o,„/7=0. It is desirable 
to select a spherical integration volume because the energy matrix may 
be nearly singular, in at least some cases, if the integration volume does 
not extend over all angles (e.g. where integration is performed over a 

20 range of angle 0 which is smaller than 0s6<tt:). 
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[0025] The formal solution for the set of coefficients a that yields 
the observed values for B and also minimizes the energy function E2 is: 

a = K-'L'iLK-'L'y'B^ QB (13) 

where Q is the backward solution matrix. Q depends only upon the 

5 locations and properties of magnetic detectors 24. The matrix LK'^L^ 
typically has a number of small eigenvalues. It is therefore generally 

necessary to regularize this matrix to obtain an approximation for Q . 

The matrix can be regularized by any suitable mathematical technique, 
such as single value decomposition. One approach to performing single 
10 value decomposition is described by Uutela, K. et al. Visualization of 
magnetoencephalographic data using minimum current estimates, 
Neurolmage v. 10, (1999) pp. 173-180. Tikhonov regularization is 
another type of procedure that may be performed to regularize LK'^LJ, 

1 5 [0026] In some alternative embodiments of the invention, is 
smaller than or equal to M (the number of sensors). In these 
embodiments the scalar potential is again approximated by the finite 
series of equation (6). The relationships between the coefficients and the 
measurements is given by: 

20 5f = I L„^a^ (14) 



These embodiments use a fitting process to identify a set of coefficients 
which provides an estimate of the scalar potential representing a best fit 
to the measured data. For example, the fitting process may identify a 
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solution (i.e. a set of coefficients a) that minimizes a weighted sum of 
errors in the fit to data B as follows: 

F2 = I (^f - ^.K.(^f - B„) (15) 

m,n-\ 

where fFis a positive definite matrix. Wrmy be chosen to be the MxM 
5 identity matrix. If so, the solution which minimizes F2 is: 

5= {l'^WLY'wB^ QB (16) 
where Q is the backward matrix for the case N<M. 

[0027] A spherical harmonic representation of the magnetic scalar 
10 potential can be used to compensate for motion of the subject's head. In 
general, any head motion can be broken down into a 3 -dimensional 
rotation and a 3-dimensional translation. Knowing the rotational and 
translational properties of the spherical harmonic components of the 
potential, one can derive a mathematical relationship which transforms 
15 the actual outputs of magnetic detectors 24 to the outputs that would 
have been measured had the subject's head remained stationary. 

[0028] Spherical harmonic functions have a simple behaviour under 

rotations. For a single value of the spherical harmonic degree /, the 

20 coefficients a^^,, are transformed by a rotation according to: 

/ 

^Imn = S^S^//^ (17) 
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where a^^^ are the transformed coefficients and i?^^ ^® elements of a 

rotation matrix. Rotation does not mix spherical harmonics with different 
values of /. R^'^ is a unitary order -(2/+1) matrix. Expressions for the 
rotation matrix are provided in D.A. Varshalovich et. al.. Quantum 
5 Theory of Angular Momentum. World Scientific, Singapore (1988) ISBN 
9971-50-107-4. 



[0029] There is no simple matrix for performing spatial translations 
of the spherical harmonic functions. Spherical harmonic functions have 
10 no special symmetry under translations. It can be shown that the shifted 
and unshifted inner harmonic coefficients are related by the matrix 
equation: 

<^A//o(") = 2 '^^M^lmQ (18) 
Im 

where the elements of f are given by: 

The elements of f can be determined numerically or by published 
methods. For example, M Danos et. al. Multiple Matrix Elements of the 
Translation Operator J. MathPhys. V.6, pp 766-778 (1965) describes 
spherical harmonic translation transformations. Certain elements of f 
20 are zero. Each element of the matrix r is a simple polynomial function 
of the components of the vector u . The polynomial coefficients may be 
precalculated and stored to facilitate fast computation of f{jA) . 



r —u 



r-u 



\r — u\ 



sin edGd(p (19) 
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[0030] To derive a matrix that transforms the outputs from 
magnetic detectors 24 to the values that those outputs would have had if 
the subject's head had not been moved, one can combine a rotation 
matrix, R , a spherical harmonic function shift matrix, f , a regularized 

5 backwards matrix Q and the forward matrix L as follows: 

5(0,0)„, « Z(0,0)„4^^^>]^f(- u)^M0.0U{R.u)^ (20) 

where B (0,0) is a vector of the corrected magnetic detector outputs; 
B {r,u)^ is a vector of the actual magnetic detector outputs for the 

subject's head at an actual position and orientation which differ from the 
1 0 reference position and orientation (0,0) by the rotation R and the 

translation u . In general, equation (20) is an approximation because Q 
is regularized in the case where N>M or coefficients for the spherical 
harmonic functions are a best fit to measured data in the case where 
N<M. 

15 

[0031] One can vary the number of terms in the spherical harmonic 
expansion used as an estimate of the scalar magnetic potential. For MEG 
in which a goal is to model electric current dipole sources at a distance of 
approximately 9 cm from magnetic detectors 24, it is desirable to use 

20 intemal spherical harmonics having values of / up to at least /=12, 

preferably /=14 and most preferably at least / =18 and at least external 
spherical harmonics having /= 2, 3 and 4. Choosing intemal spherical 
harmonic functions with /={1, 2, ...,18} and extemal spherical harmonic 
functions with /={2,3,4} results in N=381 spherical harmonic functions. 

25 Choosing intemal spherical harmonic functions with /={1, 2, ...,14} and 
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external spherical harmonic functions with /={2,3,4} results in N=309 
spherical harmonic functions. 

[0032] In general it is desirable that there be more spherical 
5 harmonic functions in the expansion of the scalar potential than there are 
magnetic detectors 24 used in determining the coefficients corresponding 
to the spherical harmonic functions. This permits a solution to be 
obtained by minimizing a total gradient as illustrated, for example, in 
equation (7). As noted above, some embodiments of the invention use 
10 fewer spherical harmonic functions than sensors. These embodiments 
can use a fitting process as illustrated, for example, by equation (15). 

[0033] The approximation of the scalar potential may include one 
or more terms in addition to spherical harmonic terms. In some 
15 embodiments the scalar potential is expressed as: 

(21) 

//=! 

where a ' is a coefficient and g{r) is a function, g (r ) may model the 
contribution to the scalar potential from one or more known extemal 
magnetic sources. For example, the magnetic field produced by a 
20 subject's heart may be modelled as a point-like magnetic source. The 
potential from a point source of strength M is given by: 

gir) = ^.3 (22) 
\r - s\ 



- 16- 

A long series of external spherical harmonic functions could be needed 
to approximate the contribution of such a source to the scalar magnetic 
potential to the same accuracy as equation (21). 

5 [0034] The processes described above for determining the 

coefficients a can also be used where the scalar potential includes one or 
more additional terms as described above. The translation and rotation 
operations would need to be modified to perform translation and rotation 
in a manner appropriate to the function g{r) . Where g{r) is a simple 

10 function as in equation (21) performing translation and/or rotations can 
be performed by way of straightforward translation and rotation 
operators. 

[0035] Certain implementations of the invention comprise computer 
1 5 processors which execute software instructions which cause the 

processors to perform a method of the invention. For example, one or 
more processors in a controller for a MEG system may implement a 
method as shown in Figure 2 by executing software instructions from a 
program memory accessible to the processor(s). The invention may also 
20 be provided in the form of a program product. The program product may 
comprise any medium which carries a set of computer-readable signals 
comprising instructions which, when executed by a computer processor, 
cause the data processor to execute a method of the invention. Program 
products according to the invention may be in any of a wide variety of 
25 forms. The program product may comprise, for example, physical media 
such as magnetic data storage media including floppy diskettes, hard disk 
drives, optical data storage media including CD ROMs, DVDs, electronic 



- 17- 

data storage media including ROMs, flash RAM, or the like or 
transmission-type media such as digital or analog communication links. 

[0036] In some embodiments of the invention, signal processing 
5 system 28 permits a user to specify how many spherical harmonic 

functions will be used in estimating the scalar potential. For example, a 
user interface to signal processing system 28 may permit a user to specify 
a maximum value of / for intemal spherical harmonic terms to be used in 
the approximation for the scalar potential and a range of values of / to be 

10 used for extemal spherical harmonic terms. For example, the user may 
choose to use extemal spherical harmonic terms for values of 1=2 to /=3 
or 1=2 to /=4 or the user may choose not to include any extemal spherical 
harmonic terms at all. The signal processing system may automatically 
determine N based upon the user input and may then select an 

15 appropriate process for determining the coefficients depending upon 
whether M<N or N<M, 

[0037] Where a component (e.g. a software module, processor, 
detector, assembly, device, circuit, etc.) is referred to above, unless 

20 otherwise indicated, reference to that component (including a reference 
to a "means") should be interpreted as including as equivalents of that 
component any component which performs the function of the described 
component (i.e., that is functionally equivalent), including components 
which are not structurally equivalent to the disclosed structure which 

25 performs the function in the illustrated exemplary embodiments of the 
invention. 



- 18- 

[0038] As will be apparent to those skilled in the art in the light of 
the foregoing disclosure, many alterations and modifications are possible 
in the practice of this invention without departing from the spirit or scope 
thereof. For example: 

• Magnetic detectors 24 may be of any suitable types. For example, 
detectors 24 could include magnetometers and/or first or second 
order magnetic gradiometers. 

• While the invention is discussed above in the context of a MEG 
system, the invention may be applied to magnetic imaging of other 
biological systems or other objects which move relative to an array 
of magnetic sensors. 

Accordingly, the scope of the invention is to be construed in accordance 
with the substance defined by the following claims. 



